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Abstract 



We study lattice QCD in the hmit that the quark mass and chemical potential 
are simultaneously made large, resulting in a controllable density of quarks which do 
not move. This is similar in spirit to the quenched approximation for zero density 
QCD. In this approximation we find that the deconfinement transition seen at zero 
density becomes a smooth crossover at any nonzero density, and that at low enough 
temperature chiral symmetry remains broken at all densities. 



^ Address after 1 Oct. 1995: Brookhaven National Laboratory, Upton, NY 11973 



1 Introduction 



Lattice QCD with a nonzero density of quarks is a difficult problem due to the fact that the 
fermion determinant is not a positive real number, and thus cannot be used as a weight for 
generating configurations by Monte Carlo methods |l|]. This is unfortunate since situations 
with a non-negligible density of quarks are interesting physically, such as the interiors of 
neutron stars or the dynamics of heavy ion collisions at RHIC. A further technical difficulty 
is that with Kogut-Susskind quarks the fermion matrix is badly ill-conditioned at nonzero 
chemical potential^, making simulations even more difficult. Hence at present, only crude 
results on very small lattices are available[^ Because of our inability to deal with the 
full problem, we may consider approximations which hopefully capture some of the essential 
features of the physics. Here we present a study of QCD at arbitrary quark density in an 
approximation where the dynamics of the quarks has been removed. This approximation is 
then analogous to the quenched approximation at zero density, an approximation which has 
provided considerable insight into the nature of QCD. 



Our idea is to simultaneously take the limits of infinite quark mass and infinite chemical 
potential in such a way that the density of quarks remains fixed at some value. This leaves 
us with quarks that can be present or absent at each lattice site, but which do not move 
in the spatial directions owing to their infinite mass. The result is a much simpler fermion 
determinant such that gauge variables can be updated to equilibrium in the background of 
a prescribed density of quarks, with little more difficulty than updating in the quenched 
approximation. 

The general idea of studying the problem in simple approximations is not new. DeGrand 
and DeTar have studied an extension of the three dimensional Potts model with an imaginary 
magnetic field, which has similar symmetry breaking as in QCD, and might be expected to lie 
in the same universality class 0. Satz has used a lowest order hopping parameter expansion 
on 8^ X 3 lattices^, which is also an approach based on very heavy quarks. 

With a chemical potential included, the lattice Dirac operator using Kogut-Susskind 
quarks is 
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where Bj contains the diagonal mass term and all the spatial hopping terms on the ith time 
slice, and Tj contains all the time direction hoppings from slice i to shce i + 

We now take the limits m — > cxo and — > 00 simultaneously. This leaves us with 2ma 
along the diagonal, and the forward hopping terms, e'^'^Tj. Each spatial point is decoupled 
from all others, and the fermion determinant is just a product of easily computed determi- 
nants on each static worldline: 



det(M) = n e'^""^"* det(P^ + CI). 



(3) 



Here is the Polyakov loop at spatial site x, Uc is the number of colors, and rit is the number 
of time slices. The coefficient of the unit matrix, C, is {2ma/e^Y\ and is the fundamental 
parameter in our approximation, through which we fix the density (We will see later that 
is the ratio of the probability that there are three quarks (in SU(3)) on a site to the 
probability that the site is empty.). 



The determinant is easily evaluated by diagonalizing P^. In SU(2), 

det(Pj + C) = + C T: + 1, 

while in SU(3) 

det(Pj + C) = + Tr + C Tr P| + 1. 



(4) 
(5) 



In SU(2) this determinant is real and positive, which reflects the fact that quarks and 
antiquarks are in equivalent representations of SU(2). Unfortunately, studying SU(2) at 
large baryon density is of limited interest since such baryons would be bosons. Neutron stars 
would be quite different in this case as would supernovae. In the realistic case of SU(3), we 
are still left with a complex determinant, albeit a much simpler one, allowing us to generate 
high statistics. 

In generating gauge configurations, the determinant in the partition function is written 
as the exponential of a sum over spatial sites 

Z = y[(iC/]ne'^~det(Pj + C)e-^« 
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We can update the spatial links with any of the standard algorithms for quenched QCD: 
Metropolis, heat bath and/or overrelaxation. The temporal links are updated with the 
Metropolis algorithm, using the magnitude of the determinant plus the gauge action as the 
weight. Thus the parts of the action involved in updating a temporal link Ut are 

S = --- + (2//)ReTrf/tS + log(|det(t/tP,^ + C)|) + ••• (7) 

where S is the sum of the "staples" and is the product of all the time direction links 
at site X except for the link being updated. As in conventional quenched QCD, successive 
Metropolis hits are easy; most of the work goes into evaluating the staples and which can 
be used unchanged in successive hits. 

Because of the phase in the determinant, for SU(3) we must estimate expectation values 
by taking the ratio 

iO) - ^ (8) 

where 6 is the phase of the determinant summed over all spatial points and ()|| indicates 
an expectation value in the ensemble of configurations weighted by the magnitude of the 
determinant. This method can be applied to the full theory; however the expectation value of 
the phase can (and typically does) become very small, so that enormous statistics are required 
to get meaningful measurements. In our case the phase does become small, however not 
prohibitively so on the lattices we study (our smallest value was {e^^)\\ = 0.02). Furthermore, 
as described above, we can produce statistics generously. 

The physical quark density is obtained from 

(„) - ' ^'"(^) 




where V is the spatial volume and /? = arit is the temporal extent of the lattice. Using 
equations ^ and ^ this becomes 
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in SU(2), and 
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in SU(3), where = Tr(P^). At C = oo the density is 0; at C = the system is saturated 
with density ric per site; C = 1 represents "half-filhng" and the density is nc/2. 

Note from Eqs. H and | that det(Pi + C) is unchanged by the replacement C 1/C, 
where in SU(3) we simuhaneously replace U —>■ U*. Thus there is a duality relation: the 
ensemble of configurations generated with coupling C is the same as that generated at 
coupling 1/C, and the density obeys p(l/C) = ric — p{C). Physically this duality reflects 
the fact that ric — l quarks on a single site behave like an antiquark on that site, so that the 
system nearly saturated with quarks, or equivalently a quark saturated system with a small 
density of holes behaving as antiquarks (at small C), behaves the same as a system with a 
small density of quarks (at large C). 

Alternatively, we may wish to know the probability that a site contains zero, one, two, or 
three (in SU(3)) quarks. Imagine separating the partition function into sectors with a fixed 
number of quarks on the particular site in question. These sectors are distinguished by their 
dependence on the chemical potential 

Z = Zo + e^^'^'Zi + e^'^'^^'Za + e^^'^^^Zg. (12) 

The fermion determinant at site x introduces (Eqs. |^ and ^ a factor of 

into the integrand of the partition function, and C contains a factor of e~^""*. Thus the 
probabilities of finding given numbers of quarks at site x are just given by the individual 
terms in this polynomial in C: 




C3 + C^Tg + CTi + 1 



+ C^Ts + CTi + l 



Note, the probability to find either zero, one, two, or three quarks is normalized to one for 
each configuration. We see from Eqs. |l^ that the ratio of the probability of three quarks 
on a site to the site being empty is just C^^. However, the value of C does not a priori 
determine the density, since the probabilities of one or two quarks on a site depend on the 
dynamics. Similarly we may calculate correlation functions {Pn{x)Pm{y)) , the probability for 
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n quarks at site x and m quarks at site such as 



^"^"^ ( (c^ + C^T^ + CT| + l) (^3 + C^Tg + CTS + l) ) ' ^^^^ 

etc. These correlations may be used to study the clustering properties of baryons in the 
model as temperature and density are varied. 

Since in our approximation the only remaining combination of chemical potential and 
quark mass is C, the heavy quark condensate <^^'\I'^ = ^^^^ is trivially related to {n). 

However we use (tpip'^ evaluated for light quarks on the generated lattices as an indicator 
of chiral symmetry breaking. This is just a probe of the nature of the gauge configurations, 
rather than a condensate of the actual quarks in the model. It represents the chiral properties 
of light valence quarks in the presence of a finite density of massive quarks. 

We can also calculate the average Polyakov loop. Not surprisingly, since the heavy quarks 
are coupled directly to the Polyakov loop, it gets a nonzero expectation value at any C other 
than zero or infinity. This does not necessarily represent deconfinement, since we have put 
in a density of quarks which can now shield the test quark represented by the Polyakov loop. 



3 Simulation Results 



Below we describe results for both SU(2) and SU(3) gauge groups; it is interesting to compare 
the two symmetries, since the nature of the phase transition is rather different in each. For 
pure SU(2) gauge theory the high temperature transition at zero density is second order [0, 
so we expect smooth behavior as the density is varied away from this point. For pure SU(3) 
(or quenched QCD), the high temperature transition is first order 0, and we expected this 
behavior to extend into the interior of the T — /i phase diagram. Surprisingly, we find this is 
apparently not the case. The first order transition appears to become a crossover at non-zero 
density, becoming quite smooth at relatively low densities. 



3.1 SU(2) 

Our initial studies were done in SU(2) where we ran at a two values of the gauge coupling, 
A/g'^ = 1.5 and 2.0 while varying the parameter C between zero and one. This actually 
corresponds to varying the density of "antiquarks" (holes relative to a quark saturated lattice) 
from zero to one per site, but as mentioned earlier, the physics is the same as for quarks. 
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Figure 1: The expectation of the density for SU(2) given by 
Eq. |10|. C = corresponds to zero density of antiquarks, and 
C = 1 corresponds to one antiquark per site. 
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Figure 2: The magnitude of the Polyakov loop for SU(2) and 
the dense static quark action given in Eq. The nonzero value 
of the Polyakov loop is due to its explicit coupling to C in the 
action and does not necessarily indicate deconfinement. 
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Figure 3: The light quark found by inverting the usual 

Kogut-Susskind fermion matrix on gauge configurations gener- 
ated with the dense static quark action given in Eq. ^. The light 
quark mass is 0.025, and we have normalized to two fiavors. 
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The second value of the gauge couphng is near the = 4 zero density temperature driven 
transition which is at roughly A/g'^ ^ 2.3, while the first corresponds to a much lower 
temperature. The simulations were done on 8^ x 4 lattices. 

In Figs. and ^ we show the density of antiquarks and the magnitude of the Polyakov 
loop as a function of the parameter C. As expected, both exhibit smooth behavior. In Fig. 
^ we show the light quark (^ipip'j (normalized to two flavors) using Kogut-Susskind quarks 
with arriq = 0.025 for various values of the density. For A/g"^ = 2.0 (near the zero density 
transition) (^ipip'^ shows a smooth decrease to a minimum that is roughly one half the zero 

density value. However, on the colder lattice there is only a slight decrease in l^^if^. Since 

the light quark is measuring the density of near zero eigenvalues of the light quark 

hopping matrix, this amounts to saying that for cold enough lattices a high density of static 
quarks does not suffice to remove the disorder in the gauge fields found at zero density. 



3.2 SU(3) 

We have run on 6^ x 2, 8^ x 2, 10^ x 2, and 6^ x 4 lattices. Typical runs include 500 
equilibration sweeps of the lattice and 4000 measuring sweeps, where in each sweep we make 
two overrelaxation updates of the spatial links and ten Metropolis updates of both spatial 
and temporal links. The average phase on the rit = 2 lattices ranges from one, at 1/C = 0.0, 
to as small as 0.04 (the 10^ x 2 lattice at 6/g^ = 4.8 and 1/C = 0.08). Since all physical 
observables are obtained from a ratio of expectation values (Eq. |^) where the numerator 
and denominator are strongly correlated, we use a jackknife procedure with ten blocks to 
estimate the errors. 

In Figs. ^ and ^ we summarize the behavior of the density and Polyakov loop magnitude 
(|P|) in the T — fi plane, on A^^^ = 2 lattices. Fig. § shows (n) as a function of Q/g'^ (temper- 
ature) and 1/C (effectively /i). In Fig. || at 1/C = 0, or zero density, we see the strong first 
order temperature induced transition at 6/g^ ~ 5.1 in \P\. As the density increases, this 
transition smooths out. 

Figs. H, 1^ and ^ show the density, \P\ and light quark (^ipip'^ along lines of constant 1/C 
on 6^ X 2, 8^ X 2 and 10^ x 2 lattices. From these plots, it appears that the first order 
transition at 1/C = becomes a smooth crossover for any nonzero value of the density. 
This is surprising, since conventional wisdom says that a nonzero discontinuity at the edge 
of a phase diagram decreases continuously to zero at some point in the interior of the phase 
diagram. However, we note that simple functions such as tanh(^^) have a discontinuity 
at h = but a crossover at any nonzero h. Since we see no systematic dependence of the 
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Figure 4: The quark density in SU(3) fornt = 2 as a function of 
6/51^ and 1/C. Note that the density is always zero at 1/C = 0. 
This plot includes results from 6^ x 2, 8^ x 2 and 10^ x 2 lattices. 
In the smoother regions of the plot some points were interpolated 
to produce a regular mesh. 
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Figure 5: The magnitude of the Polyakov loop in SU(3) as a 
function of ^/g^ and 1/C. The contour hne is where \P\ — 0.5. 
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Figure 6: The density as a function of on nt = 2 lattices. 
The curves are, from left to right, 1/C = 0.08, 1/C = 0.04 
and 1/C = 0.02. The 1/C = 0.0 curve is absent because the 
density is always zero there. The octagons are 6^ x 2 lattices, 
the squares 8^ x 2 and the diamonds 10^ x 2. The lines just 
connect the points for each lattice size. 
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Figure 7: The magnitude of the Polyakov loop averaged over 
the lattice as a function of 6/(7^ on ra^ = 2 lattices. The curves 
are, from left to right, 1/C = 0.08, 1/C = 0.04, 1/C = 0.02 and 
1/C = 0.0. The meaning of the symbols is the same as in Fig. ^ 
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Figure 8: The light quark {i^ipj as a function of 6/^"^ on = 2 
lattices. The curves are, from left to right, 1/C — 0.08, 1/C — 
0.04, 1/C = 0.02 and 1/C = 0.0. 
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Figure 9: The behavior of (e*^)|| as 1/C varies. The lattice size 
is 6^ X 4 and the gauge couphng is G/g"^ = 5.0. 



crossover on the spatial size, except for the expected decrease of the Polyakov loop magnitude 
on cold lattices, we conclude that this rounding is not a finite size effect. 

We have also made a series of runs on 6^ x 4 lattices. We began with a series of runs 
at Q/g"^ = 5.0. Since the high temperature transition at 1/C = occurs at G/g"^ = 5.1 for 
Nt = 2, this is a fairly cold lattice, with a temperature less than half that for deconfinement 
at zero density. Fig. ^ shows the behavior of the phase (e*^)||, which in this case gets as 
small as 0.02. Again, as in SU(2), we find that the light quark (tfjifj'j remains large for any 
density at this (cold) value of G/g"^, as shown in Fig. |l^. At 1/C = 1.0, the density reaches 
1.5 quarks per site, where the quarks have the largest effect on the gauge configurations. In 
Fig. |ll]we show the light quark at C = 1 as a function of Q/g'^, showing the crossover 

to the chirally symmetric phase as the temperature is raised; while smooth, the crossover is 
nonetheless quite distinct and at occurs at lower temperature than at zero density. 
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Figure 10: The light quark as a function of 1/C on 6^ x 4 

lattices, at Q/g"^ = 5.0. 



4 Remarks 



Clearly there is much more work to be done in the direction of finite density simulations. 
In this model, we have not yet addressed the question of whether these results scale in the 
continuum limit nor tried to extract physical numbers for various quantities. Perhaps the 
best hope for progress lies in the use of improved actions, with which it may be possible to 
approach the continuum limit using lattices that are small enough so that the phase problem 
is no longer hopeless. 

The mean field analysis of the Potts model in Ref. |§ showed a disappearance of the phase 
transition at large density. In the 3d Potts model the first order phase transition persists 
to some nonzero density (ie. imaginary magnetic field), with a continuously decreasing 
discontinuity in the order parameter. The hopping parameter expansion in Ref. showed 
rather smooth behavior of the Polyakov loop and (^ipip'^ on the chemical potential. 

Does the static approximation have anything to do with real QCD? Certainly the nature 
of the high temperature transition at zero density depends strongly on the presence of dy- 
namical quarks as is becoming clear from large scale simulations of full QCD[P]. However, 
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Figure 11: The light quark as a function oid/g'^ at C 

on 6"^ X 4 lattices. Recall that the chiral phase transition at 
quark density is at G/g"^ ~ 5.7 



17 



it is not a priori clear to us that a deconfinement transition or chiral symmetry restoration 
driven by high density should depend on the quarks moving, or whether the mere presence of 
the quarks would be enough. In particular, we had not expected to see the zero density first 
order transition disappear for very small quark densities, or for the signal of chiral symme- 
try restoration to vanish. This suggests that we might want to re-examine the conventional 
wisdom that a high density of quarks causes a phase transition similar to that caused by 
high temperature. 
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